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A system of equations for anisotropic hydrodynamics is derived that describes a mixture of 
anisotropic quark and gluon fluids. The consistent treatment of the zeroth, first and second mo¬ 
ments of the kinetic equations allows us to construct a new framework with more general forms of 
the anisotropic phase-space distribution functions than those used before. In this way, the main diffi- 
ciencies of the previous formulations of anisotropic hydrodynamics for mixtures have been overcome 
and the good agreement with the exact kinetic-theory results is obtained. 
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I. INTRODUCTION 

The successful description of relativistic heavy-ion col¬ 
lisions at RHIC and at the LHC in terms of relativistic 
dissipative fluid dynamics (for a recent review see [1]) has 
brought a lot of attention to studies aiming at the con¬ 
struction of the most adequate hydrodynamic framework. 
One way to achieve this task is to compare the results of 
various hydrodynamic approaches [2-12], which differ by 
the number of terms included in the formalism and by 
the values of the transport coefficients, with the results 
of the underlying microscopic kinetic theory [13-18]. The 
latter is very often used as a staring point to derive the 
specific form of the evolution equations of relativistic hy¬ 
drodynamics, however, several approximations done in 
such procedures may result in differences between the 
predictions of the kinetic theory and the hydrodynamic 
models constructed directly with its help. 

As a lot of work has been already done in this con¬ 
text for simple (i.e., one-component) fluids, the anal¬ 
ysis of mixtures has been so far quite limited [19-21], 
for some recent developments see [22-24]. One of the 
problems of the previous approaches using the concept 
of anisotropic hydrodynamics [19-21] was that they were 
based only on the zeroth and first moments of the ki¬ 
netic equations. Assuming that the distribution func¬ 
tions used in anisotropic hydrodynamics are described by 
the original Romatschke-Strickland form [25], one hnds 
underdetermined set of equations, where the number of 
unknown parameters is larger than the number of equa¬ 
tions. Consequently, to close the system of equations, in 
Refs. [19-21] the transverse-momentum scale parameters 
for quarks and gluons were taken to be equal 

In this work we develop the approach presented in 
Refs. [19-21]. We use the zeroth, first, as well as the 
second moments of the kinetic equations for quarks, anti- 


^ The transverse-momentum scale parameters can be interpreted 
also as transverse temperatures — parameters characterising 
transverse-momentum distributions. 


quarks and gluons. This allows us to use the Romatschke- 
Strickland form with more independent parameters as 
compared to the previous works. Our selection of the 
equations will be presented and discussed in greater de¬ 
tail below. Here we only mention that our equations 
include the baryon number and energy-momentum con¬ 
servation laws. This is achieved by the use of the two 
types of the Landau matching condition: the first one 
fixes the effective chemical potential /i, while the second 
one (more commonly used) fixes the effective tempera¬ 
ture T. Moreover, we use a special combination of equa¬ 
tions coming from the second moment which guarantees 
the agreement with the Israel-Stewart theory for the sys¬ 
tem approaching local equilibrium [30], see also [31-33]. 

Our approach is restricted to the one-dimensional, 
boost invariant systems [34], denoted below as (O-l-l)D 
systems. We show that the evolution equations of 
anisotropic hydrodynamics for quark and gluon fluids de¬ 
rived herein yield good agreement with the results of the 
kinetic theory. Several options for the selection of the ze¬ 
roth moment equations have been studied and the best 
choice is indicated. The new approach eliminates solu¬ 
tions with exponentially damped anisotropy parameters, 
found in [19], which do not agree with the kinetic-theory 
solutions [21]. Such solutions appeared in the cases where 
the initial conditions corresponded to oblate-prolate or 
prolate-prolate initial quark and gluon momentum dis¬ 
tribution functions. The new solutions have the same 
qualitative character for all types of the initial conditions 
(oblate-oblate, oblate-prolate, and prolate-prolate), do 
not exhibit exponential damping, and closely follow the 
kinetic-theory solutions for systems far from and close to 
local thermal equilibrium. 

The paper is written as follows: In Sec. H the kinetic 
equations for quarks, antiquarks, and gluons in the re¬ 
laxation time approximation are introduced. In Sec. HI 
the zeroth moments of the quark and antiquark equations 
are discussed in the context of the baryon number con¬ 
servation. In Sec. IV we discuss the first moment of the 
kinetic equations and analyse the energy-momentum con¬ 
servation. The second moments are discussed in Sec. V 
and our results are presented in Sec. VI. We conclude in 
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Sec. VII. Throughout the paper we use natural units and we introduce the four-vector Z defined as [11] 
the metric tensor’s signature is (-I-, —, —, —). 

2’''= 7^K,0,0, 1), 7^ = (1 - (6) 


II. KINETIC EQUATIONS 


We start our analysis with the kinetic equations for 
quarks, antiquarks and gluons written in the relaxation 
time approximation (RTA) [26-28]. They read 


P^dfj,Q^{x,p) = -p^U^ 


Q^ix,p) - Qfq{x,p) 

'^eq 


(1) 


'^eq 

where Q~^{x,p) {Q~{x,p)) is the quark (antiquark) 
phase-space distribution function, G{x,p) is the gluon 
distribution function, and Tgq is the relaxation time. The 
four-vector U describes the hydrodynamic flow in the 
Landau frame (i.e., U is defined as the eigen four-vector 
of the energy-momentum tensor). 

The quark and gluon distribution functions are as¬ 
sumed to have a generic structure [25] 


Q^(x,p) 

G{x,p) 


exp 


' ±x,-^{p-ur + ^,{p-zy 

An 


exp _ , 3 , 


where the parameters A,j and Ag define the transverse 
momentnm scale, Ag is the non-equilibrium baryon chem¬ 
ical potential of quarks, while and are the anisotropy 
parameters. In local eqnilibrium, the two anisotropy pa¬ 
rameters vanish and Eqs. (3) are reduced to the standard 
equilibrium distributions 

Qfq{x,p) = exp 
Geq(x,p) = exp > (4) 



where T is the temperature and p, is the baryon chemical 
potential. For the sake of simplicity, we assume here the 
classical Boltzmann statistics. A generalisation of the 
present results to the case of the quantum Bose-Einstein 
and Fermi-Dirac statistics is straightforward [16, 29]. 

The equilibrium distribution functions of the form (4) 
are used to define the RTA collision terms in (1) and 
(2). In this case p and T should be treated as the effec¬ 
tive baryon chemical potential and effective temperature 
that are determined by the appropriate Landau matching 
conditions. 

In addition to the flow vector U, that can be parame- 
terised in terms of the three-velocity in the standard way 
as 

=-t{l,v^,Vy,vy, 7 = (1 - (5) 


The appearance of Z is connected with the privileged 
direction of the beam axis. 

The four vectors U and Z satisfy the following normal¬ 
isation conditions 

[/2 = 1, Z2 = -1, U-Z = 0. (7) 

In the local rest frame of the fluid element, and Z>^ 
have simple forms 

C/^ = (1,0,0,0), Z'^ = (0,0,0,1). (8) 

In the (O-l-l)D case, we may further use 

= (t/r,0,0,z/T), 

= (z/T,0,0,t/T), (9) 

where r is the (longitudinal) proper time 

r = \/t‘^ — z^. (10) 


III. ZEROTH MOMENTS OF THE KINETIC 
EQUATIONS 


Integrating Eqs. (1) and (2) over three-momentum and 
inclnding the internal degrees of freedom we obtain the 
three scalar equations 

( 11 ) 

'^eq 

d^ngU^) = (12) 

'^eq 


where we have introduced the non-equilibrium and equi¬ 
librium particle densities defined by the expressions 


n„ 


\/l + Cg 

9g_ 

VI + ’ 


<eq = (13) 

ng,eq = (14) 


Note that for the (O-l-l)D system we have = d/dr 

and dgU^ = l/r. 

Instead of using Eqs. (11)-(12) we use the difference 
of the equations for quarks and antiquarks appearing in 
( 11 ), 


dr 


(n+ -Ug) + 


n+ - rig 


71~^ — 71 

^q,eq ‘'g,eq 




‘ eq 


nq ) 


(15) 
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and the following linear combination of Eqs. (11)-(12) 


dun 




dr 


dria 


Ur 


= a 


dr 


q ^9 _l_ — a) 

^eq ^eq 


(16) 


In Eq. (16) we have introduced the notation riq and n^^eq 
for the sum of the quark and antiquark densities, 


n„ = nT + n„ = 


2gq cosh(A,/A,)A3 


\/l + 


‘^9q 


^q,eq 


= + «q,eq = ^ COSh(/x/T) T . (17) 


9,eq 


It is important to note that in contrast to Eq. (15), that 
leads directly to the fundamental law of baryon number 
conservation, the use of Eq. (16) is not so well motivated. 
A straightforward treatment of Eqs. (11)-(12) suggests 
that Eq. (15) should be supplemented by the two ex¬ 
tra equations, for example, one equation for the sum of 
the quark and antiquark distributions and the other for 
the gluon distribution. It turns out, however, that the 
use of three equations obtained from the zeroth moment 
leads finally to an overdetermined system of equations. 
Therefore, we have decided to use the linear combina¬ 
tion defined by Eq. (16), where a = const is a parameter 
taken from the range 0 < a < 1. 

One may check a posteriori which value of a is the 
best by comparing the hydrodynamic results with the 
kinetic-theory results. In this way we have found that 
the best agreement is obtained for a = 1. A similar 
agreement is also obtained for the case a = 0. We note 
that the cases a = 0 and a = 1 do not introduce the cou¬ 
pling between the quark and gluon sectors at the level of 
the zeroth moment. This seems to be a desirable situa¬ 
tion, since the kinetic equations (l)-(2) treated exactly 
include the coupling between the quark and gluon sectors 
only through the energy-momentum conservation (that is 
within the first moment of Eqs. (I)-(2)) with the corre¬ 
sponding Landau matching condition, see Ref. [21]. We 
come back to the discussion of this point below Eq. (43). 


A. Baryon number conservation and the 
corresponding Landau matching 


and 




A3 

q 

■\/l + iq 


( 20 ) 


respectively. We note that the solution &(t) = boTo/r has 
the scaling form known from the Bjorken model [34]. The 
quantity bo = b{T = tq) on the right-hand side of Eq. (19) 
is the baryon number density at the initial proper time. 

Equations (19) and (20) can be solved for the chemical 
potentials Xq and g [19]. We find 


^ = sinh"^ (D) = ln[£i -f \/l + D^] (21) 


where 


DiT,Aq,^q) = 


Stt^^oTo^I -b ^q 
259 TA 3 


and 


^ = sinh ^ (— 1 = In 


D I D2 


where 


Kq{T, Aq, ^q) — 


A3 

q 


( 22 ) 


(23) 


(24) 


One should note here that the ratio D/nq = 
(Sn'^boTo)/{2gqTT^) does not depend on A^ and ^q any¬ 
more. The condition 


6eq = b, (25) 

resulting in Eq. (20), should be treated as the Landau 
matching condition that guarantees the baryon number 
conservation. It defines the effective baryon chemical po¬ 
tential /i, see Eqs. (20) and (23), in terms of the variables 
T, T, Aq and ^q. Similarly, the baryon conservation equa¬ 
tion allows us to express Xq in terms of t, T, A^ and ^q. 
In this way, in the following expressions we may com¬ 
pletely eliminate both g and Xq (the expressions depend, 
however, in the explicit way on tq and bo). 


B. Sum of the zeroth-order moments 


Equation (15) divided by a factor of three gives 
the constraint on the baryon number density 
b = (n+ — n~) /3, namely 


db ^ b 6eq ~ b 

dr T Teq 


(18) 


In order to have the baryon number conserved both the 
left- and right-hand sides of (18) should vanish. This 
yields 


&(r) 


boTo 2g, 





A3 

q 

\/l + 


(19) 


Using Eqs. (21) and (23) and the mathematical iden¬ 
tity 

cosh[sinh“^(x)] = \/l -f (26) 


we find the expressions for the quark non-equilibrium and 
equilibrium densities 

2gq yrT)^A 3 

’ 

nq,eq = 


r 


(27) 
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respectively. Using this notation we rewrite Eq. (16) in 
the form 


A. 

dr 



Vl + D^Al 

\/l + 4 


+ (1 


a) 


] 

\/l +^g ) 


+ 




Vl + E>2A3 
+ iq 


+ (1 


a) 


\/l + / 


rp3 


(a^/T-KD^y^ + (1 



(28) 


where f is the ratio of the internal degrees of freedom for 
gluons and quarks 


r = 


259 ’ 


(29) 


In the numerical calculations we use the value f = 2/3. 


IV. ENERGY-MOMENTUM CONSERVATION 
LAW 

A. Landau matching for the energy density 

The hrst moment of the sum of the kinetic equations 
(1) and (2) gives the divergence of the energy-momentum 
tensor 


Correspondingly, for the equilibrium part we find 


rpfll/ rjiflV I 

-^eq '^q',eq -^q,eQ^ 


q,eq 




( 35 ) 


where 


Tt,A = (36) 

with Ej^eq and Pi,eq being equilibrium energy density and 
pressure, respectively. The equilibrium energy densities 
for quarks and gluons are given by the expressions 


^g,eq 

^9,eq 


TT^ 

3^ 

■7r2 


yrTzu/^, 


(37) 


The energy-momentum conservation law = 0 

requires that the right-hand side of Eq. (30) vanishes. 
This is nothing else but the Landau matching condition 
for the energy-momentum conservation. This match¬ 
ing requires that the energy determined from the non¬ 
equilibrium distribution functions is the same as the en¬ 
ergy obtained with the equilibrium distribution functions 


e = e„ 


' ^eq ^9,eq T ^9,eq- 


(38) 


This leads directly to the constraint on the effective tem¬ 
perature T, 


= —U^ (TM/q + - (T/^-^ + T^n) , (30) 

"^eq 

where the energy-momentum tensor includes the 
quark and gluon contributions 


rpfll/ _ _|_ rjnfj,l/ 


(31) 


The Lorentz structure of the distribution functions for 
quarks and gluons implies the following forms of the 
energy-momentum tensors for quarks and gluons [11] 

rr = {e^ + 

-{P^,T-P^,L)Z^Z\ (32) 

where e^, Pi^L and Pi^x denote energy density, longitudi¬ 
nal and transverse pressure, respectively, and the index 
i stands for quarks {i = q) or gluons (i = 9 ). Here the 
energy densities of quarks and gluons are given by the 
expressions 


£a = 


+ 7^(4), 


TT^ 

3<7gA4 


TT^ 


n^g), 


where the function TZ{^) is defined as [10] 




2(1+0 L 


1 + 


(l-bOlan ^ a/I"' 


(33) 


(34) 


J.4 ^ A4vOT/D^7^(4)+A4f7^(e0 


1 -I- /ni + f 


B. Energy and momentum conservation 

In the (O-l-l)D case considered here the energy and 
momentum conservation takes the form 


de 

dr 


£ + Pl 


(40) 


where P^ is the sum of the longitudinal pressures for 
quarks and gluons, 


PxL = 


PnT. = 


7^i(co, 

TT^ 

^9gK 


^g,L = ^2 ^ P-Li^g), 


with TZl defined through the formula 

7^L(0 = -[2(l+0l^'(0 + l?(0]• 

This leads directly to the formula 

^ n{^,)+r 


(41) 


(42) 


(43) 


2 r 


AWi + DHi + 0)1^'(0) + ^X(i + O)I^'(O) 
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We note that Eqs. (39) and (43) couple the quark and 
gluon parameters in the similar way as these two sectors 
are coupled in the exact treatment of the kinetic equa¬ 
tions (l)-(2). Hence, the best agreement between the 
hydrodynamic equations and the kinetic theory may be 
expected if no other coupling is incorporated into the hy¬ 
drodynamic approach. This, in turn, suggests to use the 
values a = 1 or a = 0 in Eq. (16). Such a conjecture has 
been supported by our numerical calculations done with 
various values of a taken within the range 0 < a < 1. 


dent It may be taken as 



Since Eqs. (44) turned out to be very successful in the 
construction of hydrodynamic models [33], in particular, 
they are consistent with the Israel-Stewart theory for sys¬ 
tems close to equilibrium, we use this form separately for 
quarks and gluons. See also our remarks below Eq. (43). 


A. Quarks and antiquarks 


V. SECOND MOMENT OF THE KINETIC 
EQUATION 

So far we have constructed three equations, see 
Eqs. (28), (39), and (43), for five unknown functions: A,, 
Ag, ^q, ^g, and T. In order to close the system of equa¬ 
tions we need to include two extra equations. We shall 
construct them using the second moment of the kinetic 
equation. 

The second moment of the kinetic equation in the re¬ 
laxation time approximation was studied in Ref. [30], 
where a boost-invariant and cylindrically symmetric sys¬ 
tems were analyzed. In this case, it was shown that it is 
convenient to select the following equations as the basis 
for the hydrodynamic approximation 


Following the method of Ref. [31] one can derive the 
following formulas for the sum of the quark and antiquark 
distributions 

09 _ 09 _ _ ^9q^q _ /^3 

Similarly, for the equilibrium quark functions one gets 

©leq = = ©I.eq = + 


Using Eqs. (47) and (48) in (46) we find 


-^ln&i + e-29i-l-J2 In 0,7 + 0-2dj 


= ± [^_ll _iv|— [^-lll . (44) 

Teq 0/ J 3-^\Teq 0J Jj 


Here / = X,Y, Z and J = X,Y, Z are space indices. 
In the (0+l)D case the coefficients 9j have the form: 
Ox = 9 y = 0, 9z = —1/t, and 9 = Ifr. The three func¬ 
tions 0/ are obtained as contractions of the third moment 
of the distribution function with the tensor U ® I Z) I 
(where I is now the four-vector corresponding to the in¬ 
dex I). The function 0eq is obtained by the contraction 
of the third moment of the equilibrium distribution func¬ 
tion with any of the tensors U ® I ® I (the result is 
independent of I). The four-vectors U and Z are defined 
above. The four vectors X and Y in our case, where the 
transverse expansion is neglected, are given simply by the 
formulas [35] 




- ^In 




dr 1(1 + ?,)'/^ 


TeqAf^^ yr+)^ ■ 


B. Gluons 


In the case of gluons, the analogous expressions are 


09 _ 09 _ ^09 

^ ^ 1/2’ 

09 ^ 4ggA5 

7r2(I+eg)3/2’ 


A''= (0,1,0,0), y^ = (0,0,I,0). (45) 


It is important to emphasise that in our one-dimensional 
case only one out of three equations in (44) is indepen- 


^ By construction, only two equations in (44) are independent. In 
the case where the transverse flow is neglected, these two become 
degenerate. 
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0.1 0.2 0.5 1.0 2.0 5.0 10.0 



FIG. 1: (Color online) Comparison of the hydrodynamic and 
kinetic-theory results for the initial oblate-oblate configura¬ 
tions. Detailed description in the text. 


and the equilibrium functions are 


^X,eq — ®y,eq “ ^Z,eq “ ^2 ' 


Using Eqs. (50) and (51) in (46) one obtains 

= + (52) 

Equations (28), (39), and (43) together with Eqs. (49) 
and (52) represent five independent equations that allow 
us to determine five unknown functions of the proper 
time: A,, Ag, 4, 4, and T. 

VI. RESULTS 

In this Section we present our numerical results for 
four different types of initial conditions. In the first 
case the initial conditions correspond to the oblate quark 
and gluon distribution functions (oblate-oblate config¬ 
uration). In the second case the initial distribution of 
quarks is prolate, while the gluon distribution is oblate 
(prolate-oblate configuration), and in the third case the 
two distributions are prolate (prolate-prolate configura¬ 
tion). In the first three cases we set baryon number den¬ 
sity equal to zero. The effect of the non-vanishing baryon 
density is studied in the fourth case where the quark and 
gluon distributions are both oblate. 

We note that the oblate (prolate) distributions cor¬ 
respond to positive (negative) values of the anisotropy 
parameter ^ and, consequently, to the transverse pres¬ 
sure larger (smaller) than the longitudinal pressure. The 
results of the microscopic calculations indicate that the 
initial conditions in relativistic heavy-ion collisions cor¬ 
respond to the (oblate) case where Ft > Pl, as soon 
as the coherent longitudinal colour fields disappear [36- 
38]. This suggests that the oblate-oblate configuration is 
probably the most realistic one. 

In all the considered cases the initial starting (proper) 
time is r = To = 0.1 fm/c and we continue the evolution 
till r = 10 fm/c. The initial transverse-momentum pa¬ 
rameters Ag^i for quarks and gluons have been set equal 
to 1 GeV. From the Landau matching condition for the 
energy, see Eq. (39), we determine the initial tempera¬ 
ture To which is different for the cases with different ini¬ 
tial anisotropies. The value of the relaxation time used 
in this work is constant, Teq = 0.25 fm/c. A tempera¬ 
ture dependent Teq can be also used, as it is described 
in [13, 14]. However, to check the agreement with the 
kinetic theory it is enough to use a constant value. 

In the four cases presented here, the results of 
anisotropic hydrodynamics are compared with the exact 
solutions of the kinetic equations (l)-(2) that have been 
constructed for (O-l-l)D systems in Ref. [21]. We refer to 
this paper for all details connected with the exact treat¬ 
ment of Eqs. (l)-(2). Here we only emphasise that the 
initial distribution functions used in the kinetic-theory 
calculations are specified by the same initial parameters 
as those used in the hydrodynamic calculations. 
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0.1 0.2 0.5 1.0 2.0 5.0 10.0 



FIG. 2: (Color online) The same as Fig. 1 but for the prolate- FIG. 3: (Color online) The same as Fig. 1 but for the prolate- 
oblate configuration. prolate configuration. 


The results for the oblate-oblate initial configuration 
with = 1 and ^ 0,3 = 10 are shown in Fig. 1. The 
panel (a) shows the ratio of the total transverse pressure 
to the total energy density, Pt /£■ The solid line describes 
the result of the kinetic theory, while the dashed line is 
the result of anistropic hydrodynamics. The panel (b) 
shows the same ratio but for the individual quark (red) 
and gluon (black) components. The panels (c) and (d) 


show the Pl/Pt ratio for the whole system and indi¬ 
vidual components, respectively. Figures 2 and 3 show 
the same ratios as those presented in Fig. 1 but for the 
prolate-oblate (^o,<? = —0.5, ^o.g = 10) and prolate- 
prolate (^ 0,3 = —0.5, ^o,g = —0.25) initial conditions, 
respectively. The coding of the lines is the same as in 
Fig. 1. 

The results presented in Fig. 1 show that the hydro- 






















0.1 0.2 0.5 1.0 2.0 5.0 10.0 



FIG. 4: (Color online) Comparison of the results obtained 
with zero and finite initial baryon density for the oblate-oblate 
configuration. 

dynamic description agrees very well with the kinetic re¬ 
sults. One may notice, however, that the best agreement 
is achieved for total quantities that include both quarks 
and gluons. On the other hand, the worst agreement is 


obtained for gluons alone. This situation changes if we 
set a = 0 instead of a = 1, however, the overall results for 
the case a = 1 are the best. In Figs. 2 and 3 again a good 
agreement between the hydrodynamic and kinetic-theory 
results can be seen, especially for the prolate-prolate con¬ 
figurations. For the prolate-oblate case, the agreement 
is a bit worse, but we have to keep in mind that such 
an initial conhguration is extremely out of equilibrium 
— not only the two distributions are highly anisotropic 
in the momentum space, but their individual types of 
anisotropy are different. A correct description of such a 
non-equilibrium case within a hydrodynamic approach is 
challenging. 

Finally, in Fig. 4 we show the effect of finite baryon 
density. The calculations are done for the two cases: 
ba = 0 and Bq — 1 fm“^. In these two cases the ini¬ 
tial oblate-oblate configuration is assumed. In agreement 
with earlier studies we find that the effect of finite baryon 
density is very small, unless the initial baryon number 
density is extremely large {bo ^ 1 fm“^). 

We emphasise that the results shown in Figs. 1-4 have 
been obtained with a = 1 in Eq. (16), which has turned 
out to be the best choice. Similar agreement has been 
also obtained for the case a = 0. On the other hand, the 
case a = 0.5 gives much worse agreement. 


VII. CONCLUSIONS 

We have constructed a new set of equations for 
anisotropic hydrodynamics describing a mixture of 
anisotropic quark and gluon fluids. The consistent treat¬ 
ment of the zeroth, first and second moments of the ki¬ 
netic equations allows us to construct our approach with 
more general forms of the anisotropic phase-space distri¬ 
bution functions than those used in similar earlier stud¬ 
ies [19-21]. In this way, the main problems of the previ¬ 
ous formulations of anisotropic hydrodynamics for mix¬ 
tures have been overcome and the good agreement with 
the exact kinetic-theory results is obtained. 
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